rm(list=ls())
if(!is.null(sessionInfo()$otherPkgs)){
  invisible(lapply(paste0('package:', names(sessionInfo()$otherPkgs)), detach, character.only=TRUE, unload=TRUE))  
}
if (!require(pacman)){
  install.packages(pacman)
  library(pacman)
}
p_load(stargazer,
       survey,
       dplyr,
       questionr,
       srvyr,
       haven,
       jtools,
       interactions,
       ggplot2,
       coefplot,
       corrplot,
       sp,
       sf,
       ggmap,
       raster,
       geodata,
       stringdist,
       ggthemes,
       stringi,
       spdep,
       tidyr)

dta <- readRDS("../Data/clean_data.RDS")
labels <- read_dta("../Data/wave1.dta") %>% dplyr::select(starts_with("County"), Bundesland_Germ, Province_Turk) %>% unique()

########Construct labelled maps############
labshun <- attributes(labels$County_Hung)$labels
labsrom <- attributes(labels$County_Rom)$labels
labsger <- attributes(labels$Bundesland_Germ)$labels
labstur <- attributes(labels$Province_Turk)$labels


hun <- dta %>% filter(Country == 76)
rom <- dta %>% filter(Country == 142)
ger <- dta %>% filter(Country == 1358)
tur <- dta %>% filter(Country == 179)


hunadm <- st_as_sf(gadm(country='HUN', level=1,path=tempdir()))
romadm <- st_as_sf(gadm(country='ROU', level=1,path=tempdir()))
geradm <- st_as_sf(gadm(country='DEU', level=1,path=tempdir()))
turadm <- st_as_sf(gadm(country='TUR', level=1,path=tempdir()))

turnuts1 <- c("Istanbul",
              rep("West Marmara",5),
              rep("Aegean",8),
              rep("East Marmara",8),
              rep("West Anatolia",3),
              rep("Mediterranean",8),
              rep("Central Anatolia",8),
              rep("West Black Sea",10),
              rep("East Black Sea",6),
              rep("Northeast Anatolia",7),
              rep("Central East Anatolia",8),
              rep("Southeast Anatolia",9)
              )
turnuts2 <- c("Istanbul",
              rep("Tekirdag",3),
              rep("Balikesir",2),
              "Izmir",
              rep("Aydin",3),
              rep("Manisa",4),
              rep("Bursa",3),
              rep("Kocaeli",5),
              "Ankara",
              rep("Konya",2),
              rep("Antalya",3),
              rep("Adana",2),
              rep("Hatay",3),
              rep("Kirikkale",5),
              rep("Kayseri",3),
              rep("Zonguldak",3),
              rep("Kastamonu",3),
              rep("Samsun",4),
              rep("Trabzon",6),
              rep("Erzurum",3),
              rep("Agri",4),
              rep("Malatya",4),
              rep("Van",4),
              rep("Gaziantep",3),
              rep("Sanliurfa",2),
              rep("Mardin",4))
turnuts3 <-c("Istanbul",
             "Tekirdag",
             "Edirne",
             "Kirklareli",
             "Balikesir",
             "Çanakkale",
             "Izmir",
             "Aydin",
             "Denizli",
             "Mugla",
             "Manisa",
             "Afyon",
             "Kütahya",
             "Usak",
             "Bursa",
             "Eskisehir",
             "Bilecik",
             "Kocaeli",
             "Sakarya",
             "Düzce",
             "Bolu",
             "Yalova",
             "Ankara",
             "Konya",
             "Karaman",
             "Antalya",
             "Isparta",
             "Burdur",
             "Adana",
             "Mersin",
             "Hatay",
             "K. Maras",
             "Osmaniye",
             "Kinkkale",
             "Aksaray",
             "Nigde",
             "Nevsehir",
             "Kirsehir",
             "Kayseri",
             "Sivas",
             "Yozgat",
             "Zinguldak",
             "Karabük",
             "Bartın",
             "Kastamonu",
             "Çankiri",
             "Sinop",
             "Samsun",
             "Tokat",
             "Çorum",
             "Amasya",
             "Trabzon",
             "Ordu",
             "Giresun",
             "Rize",
             "Artvin",
             "Gümüshane",
             "Erzurum",
             "Erzincan",
             "Bayburt",
             "Agri",
             "Kars",
             "Iğdır",
             "Ardahan",
             "Malatya",
             "Elazığ",
             "Bingöl",
             "Tunceli",
             "Van",
             "Mus",
             "Bitlis",
             "Hakkari",
             "Gaziantep",
             "Adiyaman",
             "Kilis",
             "Sanliurfa",
             "Diyarbakir",
             "Mardin",
             "Batman",
             "Sirnak",
             "Siirt")
turnuts <- data.frame(nuts1 = turnuts1, nuts2 = turnuts2, nuts3=turnuts3)

#Create nuts2 polygons of Turkey
subregions <- unique(turnuts2)
subregionssf <- list()
for(i in 1:length(subregions)){
  temp <- turnuts %>% filter(nuts2 == subregions[i])
  subregionssf[[i]] <- turadm[turadm$NAME_1 %in% temp$nuts3, ] %>% 
    st_buffer(0.5) %>% 
    st_union() %>% 
    st_sf() %>% 
    mutate(NAME_2 = subregions[i])  
}
turadm2 <- do.call(rbind,subregionssf)
regions <- unique(turnuts1)
regionssf <- list()
for(i in 1:length(regions)){
  temp <- turnuts %>% filter(nuts1 == regions[i])
  regionssf[[i]] <- turadm[turadm$NAME_1 %in% temp$nuts3, ] %>% 
    st_buffer(0.5) %>% 
    st_union() %>% 
    st_sf() %>% 
    mutate(nuts1 = regions[i])
}
turadm1 <- do.call(rbind,regionssf)
tur <- tur %>%
  mutate(Province_Turk = recode(Province_Turk,!!!as.list(setNames(names(labstur),labstur))))
turnuts$nuts3 <- stri_trans_general(turnuts$nuts3,"Latin-ASCII")
tur$Province_Turk <- stri_trans_general(tur$Province_Turk,"Latin-ASCII")
tur$Province_Turk[tur$Province_Turk=="Afyonkarahisar"] <- "Afyon"
tur$Province_Turk[tur$Province_Turk=="Kahramanmaras"] <- "K. Maras"
turnuts$nuts3[turnuts$nuts3=="Zinguldak"] <- "Zonguldak"
turnuts$nuts3[turnuts$nuts3=="Kinkkale"] <- "Kirikkale"
turmatch <- amatch(tur$Province_Turk,turnuts$nuts3,method="lcs",maxDist=3)
tur$matchid <- turmatch
turnuts$matchid <- seq_along(turnuts3)
tur <- merge(tur,turnuts,by="matchid",all.x=T)
#print(tur %>% group_by(nuts1) %>%summarise(n=n()),n=26)

###Average Concern by Region for Each Country###
hunavgcon <- hun %>% group_by(County_Hung) %>% summarise(avgcon = mean(territ_loss_concern,na.rm=T),romaniatransylvania=mean(romaniatransylvania),
                                                         upperhungary = mean(upperhungary),
                                                         vojevodina = mean(vojevodina),
                                                         transcarpathia = mean(transcarpathia),n=n()) %>% 
  mutate(County_Hung = recode(County_Hung,!!!as.list(setNames(names(labshun),labshun))))


romavgcon <- rom %>% group_by(County_Rom) %>% summarise(avgcon = mean(territ_loss_concern,na.rm=T),cadrilaterul=mean(cadrilaterul),
                                                        bucovina = mean(bucovina),
                                                        bessarabia = mean(bessarabia),n=n()) %>% 
  mutate(County_Rom = recode(County_Rom,!!!as.list(setNames(names(labsrom),labsrom))),
         most_named = case_when(
           bessarabia > bucovina & bessarabia > cadrilaterul ~ "Bessarabia",
           bucovina > bessarabia & bucovina > cadrilaterul ~ "Bucovina",
           cadrilaterul > bessarabia & cadrilaterul > bucovina ~ "Cadrilaterul",
           TRUE ~ "Other"
         ))
geravgcon <- ger %>% group_by(Bundesland_Germ) %>% summarise(avgcon = mean(territ_loss_concern,na.rm=T),poland=mean(poland),
                                                             france = mean(france),
                                                             russia = mean(russia),
                                                             sudetenland = mean(sudetenland),n=n()) %>% 
  mutate(Bundesland_Germ = recode(Bundesland_Germ,!!!as.list(setNames(names(labsger),labsger))))


turavgcon <- tur %>% group_by(nuts1) %>% summarise(avgcon = mean(territ_loss_concern,na.rm=T),greece=mean(greece),
                                                           regimesupport = mean(like_akp,na.rm=T),
                                                           iraq = mean(iraq),
                                                           syria = mean(syria),
                                                           bulgaria = mean(bulgaria),
                                                           otherbalkans = mean(otherbalkans),n=n()) %>% 
  rowwise() %>%
  mutate(most_named = names(.)[3:7][which.max(c(greece,iraq,syria,bulgaria,otherbalkans))],most_named = recode(most_named,greece="Greece",iraq="Iraq",syria="Syria",bulgaria="Bulgaria",otherbalkans="Balkans"))
  


###Merge with admin shape files###

hunadm$NAME_1 <- stri_trans_general(hunadm$NAME_1,"Latin-ASCII")
hunavgcon$County_Hung <- stri_trans_general(hunavgcon$County_Hung,"Latin-ASCII")
hunadm$soundex <- phonetic(hunadm$NAME_1,method = c("soundex"), 
                           useBytes = FALSE)
hunavgcon$soundex <- phonetic(hunavgcon$County_Hung,method = c("soundex"), 
                           useBytes = FALSE)
hungeo <- merge(hunadm,hunavgcon,by="soundex")

romadm$NAME_1 <- stri_trans_general(romadm$NAME_1,"Latin-ASCII")
romavgcon$County_Rom <- stri_trans_general(romavgcon$County_Rom,"Latin-ASCII")
romadm$soundex <- phonetic(romadm$NAME_1,method = c("soundex"), 
                           useBytes = FALSE)
romavgcon$soundex <- phonetic(romavgcon$County_Rom,method = c("soundex"), 
                              useBytes = FALSE)
romgeo <- merge(romadm,romavgcon,by="soundex")

geradm$NAME_1 <- stri_trans_general(geradm$NAME_1,"Latin-ASCII")
geravgcon$Bundesland_Germ <- stri_trans_general(geravgcon$Bundesland_Germ,"Latin-ASCII")
geradm$soundex <- phonetic(geradm$NAME_1,method = c("soundex"), 
                           useBytes = FALSE)
geravgcon$soundex <- phonetic(geravgcon$Bundesland_Germ,method = c("soundex"), 
                              useBytes = FALSE)

gergeo <- merge(geradm,geravgcon,by="soundex")
turgeo <- merge(turadm1,turavgcon,by="nuts1")
################Figure 2###########################
hungeo %>% mutate(`Average Concern` = avgcon) %>% ggplot(.,aes(fill=`Average Concern`,lab=NAME_1)) + geom_sf() + theme_map() + scale_fill_distiller(direction = 1) + theme(legend.position = "right")
ggsave("../Figures/hunavgcon.png")
romgeo %>% mutate(`Average Concern` = avgcon) %>% ggplot(.,aes(fill=`Average Concern`,lab=NAME_1)) + geom_sf() + theme_map() + scale_fill_distiller(direction = 1) + theme(legend.position = "right")
ggsave("../Figures/romavgcon.png")
gergeo %>% mutate(`Average Concern` = avgcon) %>% ggplot(.,aes(fill=`Average Concern`,lab=NAME_1)) + geom_sf() + theme_map() + scale_fill_distiller(direction = 1) + theme(legend.position = "right")
ggsave("../Figures/geravgcon.png")
turgeo %>% mutate(`Average Concern` = avgcon) %>% ggplot(.,aes(fill=`Average Concern`,lab=nuts1)) + geom_sf() + theme_map() + scale_fill_distiller(direction = 1) + theme(legend.position = "right")
ggsave("../Figures/turavgcon.png")
################Figure 3###########################
hungeo %>% mutate(`Transylvania (Erdely)` = romaniatransylvania) %>% ggplot(.,aes(fill=`Transylvania (Erdely)`,lab=NAME_1)) + geom_sf() + theme_map() + scale_fill_distiller(direction = 1) + theme(legend.position = "right")
ggsave("../Figures/hungeo.png")
romgeo %>% mutate(`Moldova (Bessarabia)` = bessarabia) %>% ggplot(.,aes(fill=`Moldova (Bessarabia)`,lab=NAME_1)) + geom_sf() + theme_map() + scale_fill_distiller(direction = 1) + theme(legend.position = "right")
ggsave("../Figures/romgeo.png")
gergeo %>% mutate(`Poland (Silesia)` = poland) %>% ggplot(.,aes(fill=`Poland (Silesia)`,lab=NAME_1)) + geom_sf() + theme_map() + scale_fill_distiller(direction = 1) + theme(legend.position = "right")
ggsave("../Figures/gergeo.png")
turgeo %>% mutate(`Greece (Aegean Islands)` = greece) %>% ggplot(.,aes(fill=`Greece (Aegean Islands)`,lab=nuts1)) + geom_sf() + theme_map() + scale_fill_distiller(direction = 1) + theme(legend.position = "right")
ggsave("../Figures/turgeo.png")

